Frequency-domain seismic modeling with the MUMPS

Seismic imaging has many applications in civil engineering, risk hazard, waste storage monitoring, oil exploration and tectonophysics. Among the different seismic imaging methods the full waveform inversion aims to exploit the full information content of seismic data through the complete resolution of the partial differential wave equation. Once the wave equation has been solved for each seismic source of the experiment, the imaging process consists of estimating the subsurface properties embedded in the coefficients of the wave equation from their solutions by solving a local optimization problem (gradientbased method). Solving the wave equation for multiple right-hand sides (i.e., seismic sources) is by far the most computationally expensive task in the seismic imaging process, and hence a careful attention should be paid to the numerical scheme that is used to solve the wave equation. The wave equation can be essentially written in the time-space or in the frequency-space domains although other alternatives are possible. In the latter case, solving the wave equation describes an elliptic boundary-value problem which requires the resolution of a large and sparse system of linear equations. One advantage of the frequency domain approach is that, if a LU factorization of the sparse matrix can be performed, solutions for multiple right hand sides can be obtained efficiently, a key issue in seismic imaging which involves several thousands of seismic sources. We shall present some results obtained with MUMPS to perform seismic modeling in the frequency domain, for simulations involving several millions of unknowns. Preliminary results obtained with the block low-rank version of MUMPS will highlight the decrease of the operation count, the memory demand and volume of communications, which can be achieved by taking advantage of the low-rank property of the elliptic wave-equation operator.